home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / math / alged34.zip / ALGEDSRC.ZIP / ALGRAPH.C < prev    next >
C/C++ Source or Header  |  1996-06-06  |  26KB  |  839 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. #include <graphics.h>
  11. #include <dos.h>
  12. /*
  13.   xs is size
  14.   xo is origin (left edge)
  15.   xz is maxx/xs (number of pixels per unit length)
  16.   np is number of points
  17. */
  18. static node* defx;     // default fx
  19. static double xs=0,xo,ys,yo,xz,yz,np,ts,to;
  20. static double err,step;
  21. static node ivar[10];       /* independent variables */
  22. static int gm,polar,nv,rain,fill;
  23. static int co1,co2,co3;    /* colors */
  24. extern char gdriver[80];
  25. extern int gmode,psz,pst;
  26. /*-----------------------------------------------------------------
  27.   Palette variables
  28. */
  29. typedef struct { unsigned char r,g,b; } palen;
  30. static palen pal[480];
  31. static int v256;
  32. /*-----------------------------------------------------------------
  33.   3d variables
  34. */
  35. static int in3d,vix=1;
  36. static double focalen = 5, cont;
  37. static double cx,cy,cz;    // camera location
  38. static double ax,ay,az;    // camera orientation
  39. static double m[4][4];     // transform matrix
  40. extern double lightx,lighty,lightz;  // see algfile.c
  41. #define EPS 1E-200
  42. /* These are the zooming ratios.... rr1 is (1-rr)/2 */
  43. #define rr 0.7
  44. #define rr1 0.15
  45. #define tt (M_PI/72)
  46.  
  47. typedef struct {       // save previous row data for polyfill
  48.   short x2,y2;          // screen location
  49. } point2d;
  50.  
  51. typedef struct {       // save previous row data for polyfill
  52.   double x,y,z;         // 3d location
  53. } point3d;
  54.  
  55. /*-----------------------------------------------------------------
  56.    find var
  57. */
  58. double getvar(char *name) { int i;
  59.   for (i=0; i<nv; ++i)
  60.     if (!strcmp(name,ivar[i].name)) return ivar[i].value;
  61.   return 0.0;
  62. }
  63.  
  64. /*-----------------------------------------------------------------
  65.   min and max
  66.   */
  67. double mind(double x, double y) {
  68.   return min(x,y);
  69. }
  70. double maxd(double x, double y) {
  71.   return max(x,y);
  72. }
  73.  
  74. /*-----------------------------------------------------------------
  75.    gray_palette
  76. */
  77. void gray_palette()
  78. {
  79.   int i;
  80.   for (i=0; i<120; ++i)   pal[i].r = pal[i].g = pal[i].b = i/2+3;
  81.   for (i=120; i<240; ++i) pal[i].r = pal[i].g = pal[i].b = 123-i/2;
  82.   for (i=0; i<240; ++i)   pal[i+240] = pal[i];
  83. }
  84.  
  85. /*-----------------------------------------------------------------
  86.    wave_palette
  87. */
  88. #define cwav(r) 31*(sin((i+r)*M_PI/120)+1)
  89. // #define cwav(r) max(min(abs((i+r)%240-120)-40,40),0)*63/40
  90. void wave_palette()
  91. {
  92.   int i;
  93.   for (i=0; i<240; ++i) {
  94.     pal[i].r = cwav(0);
  95.     pal[i].g = cwav(80);
  96.     pal[i].b = cwav(160);
  97.   }
  98.   for (i=0; i<240; ++i)   pal[i+240] = pal[i];
  99. }
  100.  
  101. /*  positive modulo  */
  102. long pmod(long x,long y) {
  103.   if (y==0) return x;
  104.   if (y<0) y=-y;
  105.   x = x % y;
  106.   if (x < 0) x += y;
  107.   return x;
  108. }
  109.  
  110. /*-----------------------------------------------------------------
  111.    set_palette - set the top 240 colors of the vga palette
  112. */
  113. void set_palette(int i)
  114. {
  115.  struct REGPACK r;
  116.  i = pmod(i,240);
  117.  r.r_ax = 0x1012;
  118.  r.r_bx = 0x0010;     /* start */
  119.  r.r_cx = 0x00F0;     /* length */
  120.  r.r_es = FP_SEG(pal);
  121.  r.r_dx = FP_OFF(pal) + i*sizeof*pal;
  122.  intr(0x10, &r);
  123. }
  124.  
  125. /*-----------------------------------------------------------------
  126.    eval - return value of p, using ivar
  127. */
  128. double eval(node *p) {
  129.   double u,v;
  130.   if (!p) return 0;
  131.   switch (p->kind) {
  132.   case VAR: return getvar(p->name);
  133.   case NUM: return p->value;
  134.   case ADD: return eval(p->lf) + eval(p->rt);
  135.   case SUB: return eval(p->lf) - eval(p->rt);
  136.   case MUL: return eval(p->lf) * eval(p->rt);
  137.   case DIV:
  138.     v = eval(p->rt);
  139.     if (fabs(v)<err) return 0;
  140.     return eval(p->lf) / v;
  141.   case EXP:
  142.     u = eval(p->lf);
  143.     v = eval(p->rt);
  144.     if (u<err && !whole(v)) modf(v,&v);    // result would be imaginary
  145.     if (fabs(v * log10(fabs(u)+1)) > 300) return 0;
  146.     if (!v) return 0;
  147.     return pow(u,v);
  148.   case EQU: return mind(eval(p->lf),eval(p->rt));
  149.   case FUN:
  150.     if (!strcmp(p->name,"sin")) return sin(eval(p->lf));
  151.     if (!strcmp(p->name,"cos")) return cos(eval(p->lf));
  152.     if (!strcmp(p->name,"tan")) return tan(eval(p->lf));
  153.     if (!strcmp(p->name,"asin")) {
  154.       v = eval(p->lf);
  155.       if (v>1 || v<-1) return 0;
  156.       return asin(v);
  157.     }
  158.     if (!strcmp(p->name,"acos")) {
  159.       v = eval(p->lf);
  160.       if (v>1 || v<-1) return 0;
  161.       return acos(v);
  162.     }
  163.     if (!strcmp(p->name,"atan")) return atan(eval(p->lf));
  164.     if (!strcmp(p->name,"sinh")) return sinh(eval(p->lf));
  165.     if (!strcmp(p->name,"cosh")) return cosh(eval(p->lf));
  166.     if (!strcmp(p->name,"tanh")) return tanh(eval(p->lf));
  167.     if (!strcmp(p->name,"rand")) return eval(p->lf)*rand()/RAND_MAX;
  168.     if (!strcmp(p->name,"sign")) {
  169.       v = eval(p->lf);
  170.       return v<0?-1:v>0?1:0;
  171.     }
  172.     if (!strcmp(p->name,"ln")) {
  173.       v = eval(p->lf);
  174.       if (v<=0) return 0;
  175.       return log(v);
  176.     }
  177.     if (!strcmp(p->name,"log")) {
  178.       v = eval(p->lf);
  179.       if (v<=0) return 0;
  180.       return log10(v);
  181.     }
  182.     if (!strcmp(p->name,"abs")) return fabs(eval(p->lf));
  183.     if (!strcmp(p->name,"r") && p->nump==2) return hypot(eval(p->lf),eval(p->rt));
  184.     if (!strcmp(p->name,"min") && p->nump==2) return mind(eval(p->lf),eval(p->rt));
  185.     if (!strcmp(p->name,"max") && p->nump==2) return maxd(eval(p->lf),eval(p->rt));
  186.     if (!strcmp(p->name,"atan2") && p->nump==2) {
  187.       u = eval(p->lf);
  188.       v = eval(p->rt);
  189.       if (v==0 && u==0) return 0;
  190.       return atan2(u,v);
  191.     }
  192.     if (!strcmp(p->name,"mod") && p->nump==2) return fmod(eval(p->lf),eval(p->rt));
  193.   default: return 0;
  194.   }
  195. }
  196.  
  197. /*-----------------------------------------------------------------
  198.    reset window,  use a weird number for xs to get better safety
  199. */
  200. void resetx() { int i;
  201.   xs = 7-M_PI/13; ys = xs*0.75;
  202.   xo=xs/-2; yo=ys/-2; np=25;
  203.   err = pow(10,-sigdig);  step=1.0;
  204.   for (i=0; i<10; ++i) ivar[i].value = 0;
  205.   cz=2.5; cy=2.58819; cx=9.33013;   // camera location
  206.   ax=.261799; ay=M_PI/2+ax; az=0;   // camera orientation
  207.   ts=xs; to=xo;      vix=1;
  208.   co1=15; co2=2; co3=1; rain=0;
  209.   focalen=5; fill=0; polar=0; cont=0.2;
  210. }
  211. /*-----------------------------------------------------------------
  212.    add a new var
  213. */
  214. void addvar(char *name) { int i;
  215.   for (i=0; i<nv; ++i) if (!strcmp(name,ivar[i].name)) break;
  216.   if (i==nv && nv<10) strcpy(ivar[nv++].name,name);
  217. }
  218. /*-----------------------------------------------------------------
  219.    look for the first var in p
  220. */
  221. void setivar(node *p) { int i;
  222.   if (p->kind==VAR) addvar(p->name);
  223.   for (i=0; i<p->nump; ++i) setivar(p->parm[i]);
  224. }
  225.  
  226. /*-----------------------------------------------------------------
  227.    show help file
  228. */
  229. void showghelp() {
  230.   FILE *f;
  231.   int i,c=0;
  232.   static char s[85];
  233.  
  234.   strcpy(s,"alged");
  235.   strcat(s,lang);
  236.   strcat(s,".hlq");
  237.   f = fopen(s,"r");
  238.   if (!f) { printf(msg[16],s);
  239.     pause; return;
  240.   }
  241.   restorecrtmode();
  242.   if (ti.screenheight>25) textmode(64);    /* ega 43 line mode */
  243.   textattr(norm);
  244.   clrscr();
  245.   i = ti.screenheight-1;
  246.   while (!feof(f)) {
  247.     printf(fgets(s,80,f));
  248.     if (!--i) {
  249.       i = ti.screenheight-4;
  250.       while (!(c=getch()));
  251.       if (c==27) break;
  252.     }
  253.   }
  254.   if (c!=27) getch();
  255.   fclose(f);
  256.   setgraphmode(gmode);
  257. }
  258. /*-----------------------------------------------------------------
  259.    rotate vars     starting at index i
  260. */
  261. void rot_var(int i) {
  262.   node t;
  263.   if (nv < 2) return;
  264.   t = ivar[i];
  265.   t.value = 0;
  266.   for (++i; i<nv; ++i) {
  267.     ivar[i-1] = ivar[i];
  268.     ivar[i-1].value = 0;
  269.   }
  270.   ivar[nv-1] = t;
  271.   strcpy(defx->name,ivar[0].name);    // default fx is ivar[0]
  272. }
  273. /*-----------------------------------------------------------------
  274.    handle key
  275. */
  276. void handlekey(char c) {
  277.   switch (c) {
  278.     case ';': showghelp(); break;     /* F1 */
  279.     case 72: yo+=ys/4; break;
  280.     case 75: xo-=xs/4; break;
  281.     case 77: xo+=xs/4; break;
  282.     case 80: yo-=ys/4; break;
  283.     case 82: xo+=xs*rr1; yo+=ys*rr1; xs*=rr; ys*=rr; break;
  284.     case 83: xs/=rr; ys/=rr; xo-=xs*rr1; yo-=ys*rr1; break;
  285.     case 71: np/=rr;  break;
  286.     case 79: if (np>6) np*=rr;  break;
  287.     case 73: ys/=rr; yo-=ys*rr1; break;
  288.     case 81: yo+=ys*rr1; ys*=rr; break;
  289.     case 'd': resetx(); break;
  290.     case 'a': polar = ++polar%3;
  291.               if (!in3d && polar>1) polar=0; break;
  292.     case 'q': ivar[1].value -= step; break;
  293.     case 'w': ivar[2].value -= step; break;
  294.     case 'e': ivar[3].value -= step; break;
  295.     case 'r': ivar[4].value -= step; break;
  296.     case 't': ivar[5].value -= step; break;
  297.     case 'y': ivar[6].value -= step; break;
  298.     case 'u': ivar[7].value -= step; break;
  299.     case 'i': ivar[8].value -= step; break;
  300.     case 'o': ivar[9].value -= step; break;
  301.     case 'p': step *= sqrt(.1);      break;
  302.     case '1': ivar[1].value += step; break;
  303.     case '2': ivar[2].value += step; break;
  304.     case '3': ivar[3].value += step; break;
  305.     case '4': ivar[4].value += step; break;
  306.     case '5': ivar[5].value += step; break;
  307.     case '6': ivar[6].value += step; break;
  308.     case '7': ivar[7].value += step; break;
  309.     case '8': ivar[8].value += step; break;
  310.     case '9': ivar[9].value += step; break;
  311.     case '0': step /= sqrt(.1);      break;
  312.     case 'V':
  313.       if (in3d) rot_var(1);
  314.       if (++vix >= nv) vix = 1;  // variable index counter
  315.       if (!in3d || vix==1) rot_var(0);
  316.       break;
  317.   }
  318. }
  319.  
  320. /*-----------------------------------------------------------------
  321.    MAIN GRAPH ROUTINE
  322.  
  323.    fx is an optional scaling function for the x coordinate, and fy
  324.    is a required function for the y coord.  In case of polar coordinates
  325.    the x is theta and the y is radius.
  326. */
  327. #define rline(x,y,a,b) line(xz*(x-xo),yz*(ys+yo-(y)),xz*(a-xo),yz*(ys+yo-(b)))
  328.  
  329. void graph(node *fx,node *fy) {
  330.   int gd=DETECT,c,i;
  331.   double t,x1,y1,x2,y2,z;
  332.   char ss[60];
  333.  
  334.   nv=0;
  335.   strcpy(ivar[1].name,"y");
  336.   setivar(fy);
  337.   if (fx) setivar(fx);
  338.   if (!nv) addvar("x");
  339.   if (!defx) defx = newvar("x");
  340.   if (!fx) {
  341.     fx = defx;
  342.     strcpy(defx->name,ivar[0].name);
  343.   }
  344.   if (!xs) {
  345.     if (*gdriver && *gdriver!='*') {    // customized bgi file!
  346.       gd=installuserdriver(gdriver,NULL);
  347.       gm=gmode;
  348.       if (!!(i=graphresult())) {
  349.         printf(msg[17],i,grapherrormsg(i));
  350.         pause;
  351.         return;
  352.       }
  353.     }
  354.     initgraph(&gd,&gm,"");
  355.     if (!!(i=graphresult())) {
  356.       printf(msg[17],i,grapherrormsg(i));
  357.       pause;
  358.       return;
  359.     }
  360.     if (!*gdriver) gmode = gm;         // auto detect and auto mode
  361.     if (gm != gmode) setgraphmode(gmode);
  362.     resetx();
  363.     if (psz<2) { psz = getpalettesize()-1; pst=1; }
  364.     v256 = (psz+pst>17 || getpalettesize()>17);
  365.   }
  366.   else setgraphmode(gmode);
  367.   while (1) {
  368.     xz = getmaxx()/xs;
  369.     yz = getmaxy()/ys;
  370.  
  371.     if (in3d && graph3d(fx,fy)) break;
  372.     cleardevice();
  373.  
  374. //    setlinestyle(DOTTED_LINE,0,1);
  375.     setcolor(co2);
  376.     rline(xo,0,xo+xs,0);
  377.     rline(0,yo,0,yo+ys);
  378.     rline(.1,1,-.1,1);
  379.     rline(.1,-1,-.1,-1);
  380.     rline(1,.1,1,-.1);
  381.     rline(-1,.1,-1,-.1);
  382.     rline(.1,2,-.1,2);
  383.     rline(.1,-2,-.1,-2);
  384.     rline(2,.1,2,-.1);
  385.     rline(-2,.1,-2,-.1);
  386.  
  387.     rline(1.1,1,.9,1);
  388.     rline(1,1.1,1,.9);
  389.     rline(-1.1,1,-.9,1);
  390.     rline(-1,1.1,-1,.9);
  391.     rline(1.1,-1,.9,-1);
  392.     rline(1,-1.1,1,-.9);
  393.     rline(-1.1,-1,-.9,-1);
  394.     rline(-1,-1.1,-1,-.9);
  395.  
  396. //    setlinestyle(SOLID_LINE,0,1);
  397.     setcolor(co1);
  398.     sprintf(ss,"%g  ",yo+ys);  outtextxy(getmaxx()/2,0,ss);
  399.     sprintf(ss,"%g  ",yo   );  outtextxy(getmaxx()/2,getmaxy()-8,ss);
  400.     sprintf(ss,"%g  ",xo+xs);  outtextxy(getmaxx()-50,getmaxy()/2,ss);
  401.     sprintf(ss,"%g  ",xo   );  outtextxy(0,getmaxy()/2,ss);
  402.     sprintf(ss,"%s %g  ",msg[18],step);
  403.     if (polar) strcat(ss,msg[19]);
  404.     strcat(ss,msg[36]);
  405.     outtextxy(0,0,ss);
  406.     sprintf(ss,"%s =[%8.4f %8.4f]  ",ivar[0].name,xo,xo+xs);
  407.     outtextxy(0,10,ss);
  408.     for (i=1; i<nv; ) {
  409.       sprintf(ss,"%s = %8.4f  ",ivar[i].name,ivar[i].value);
  410.       outtextxy(0,++i*10,ss);
  411.     }
  412.     srand(2);
  413.     for (t=xo; t<=xs+xo; t+=xs/np) {
  414.       ivar[0].value = t;
  415.       if (polar) ivar[0].value = (2*(t-xo)/xs-1)*M_PI;
  416.       x2 = eval(fx);
  417.       y2 = eval(fy);
  418.       if (polar) {
  419.         z = x2;
  420.         x2 = y2*cos(z);
  421.         y2 = y2*sin(z);
  422.       }
  423.       x2 = xz*(x2-xo);
  424.       y2 = yz*(ys+yo-y2);
  425.       if (x2>30000) x2=30000;
  426.       if (x2<-30000) x2=-30000;
  427.       if (y2>30000) y2=30000;
  428.       if (y2<-30000) y2=-30000;
  429.       if (t!=xo) line(x1,y1,x2,y2);
  430.       x1=x2; y1=y2;
  431.       if (kbhit()) break;
  432.     }
  433.     c = getch();
  434.     while (kbhit()) c = getch();
  435.     if (c==27) break;                      /* escape */
  436.     in3d=1;
  437.     if (c=='g' && graph3d(fx,fy)) break;   /* 3d graphics */
  438.     in3d=0;
  439.     handlekey(c);
  440.   }
  441.   restorecrtmode();
  442.   if (ti.screenheight>25) textmode(64);    /* ega 43 line mode */
  443. }
  444.  
  445. /*-----------------------------------------------------------------
  446.    compute matrix
  447. */
  448. static void computematrix() {
  449.   double sa,sb,sc,ca,cb,cc; int i;
  450.  
  451.   //  Precompute some sines and cosines
  452.   ca=cos(ax); sa=sin(ax);
  453.   cb=cos(ay); sb=sin(ay);
  454.   cc=cos(az); sc=sin(az);
  455.  
  456.   //  Compute rotation transformation
  457.   m[0][0]=cb*cc-sb*sa*sc;   /* R = Rb.Ra.Rc */
  458.   m[0][1]=ca*sc;
  459.   m[0][2]=sb*cc+cb*sa*sc;
  460.   m[1][0]=-cb*sc-sb*sa*cc;
  461.   m[1][1]=ca*cc;
  462.   m[1][2]=-sb*sc+cb*sa*cc;
  463.   m[2][0]=-sb*ca;
  464.   m[2][1]=-sa;
  465.   m[2][2]=cb*ca;
  466.  
  467.   //  Compute translation transformation
  468.   for (i=0; i<3; ++i) m[i][3]=-cx*m[i][0]-cy*m[i][1]-cz*m[i][2];
  469.  
  470.   for (i=0; i<4; ++i) m[3][i]=m[2][i]/focalen;
  471. }
  472.  
  473. static void view3d(double dx,double dy,double dz,double *xx,double *yy) {
  474.   double x,y,z,t;
  475.   x = dx*m[0][0] + dy*m[0][1] + dz*m[0][2] + m[0][3];
  476.   y = dx*m[1][0] + dy*m[1][1] + dz*m[1][2] + m[1][3];
  477.   z = dx*m[3][0] + dy*m[3][1] + dz*m[3][2] + m[3][3];
  478.  
  479.   //  Check if the point is behind the visual cone.
  480.   t = hypot(x,y)/5;    // 5 is slope of the cone
  481.   if (z < t) z=t;      // push point to the cone
  482.   if (z < EPS) z=EPS;  // don't divide by zero
  483.   *xx = x/z;
  484.   *yy = y/z;
  485. }
  486.  
  487. /*-----------------------------------------------------------------
  488.    move camera forward or backward
  489. */
  490. void movecam(double spd) {
  491.   double y,d;
  492.   y=sin(ax);
  493.   d=cos(ax);
  494.   cz+=spd*cos(ay)*d;
  495.   cy-=spd*y;
  496.   cx-=spd*sin(ay)*d;
  497. }
  498. /*-----------------------------------------------------------------
  499.    rotate camera - with constant radius v=-1,1 vertically or 3,5 horiz
  500. */
  501. void spincam(int v) {
  502.   double r,t,a;
  503.   r = hypot(cx,cz);
  504.   if (v<2) {              // vertical
  505.     t = atan2(cy,r) + tt*v;
  506.     if (fabs(t)>M_PI/2) return;
  507.     a = hypot(r,cy);
  508.     cy = sin(t)*a;
  509.     a *= cos(t)/r;
  510.     cx *= a; cz *= a;
  511.     ax+=tt*v;
  512.   }
  513.   else {               // horizontal
  514.     t = atan2(cz,cx);
  515.     t+=tt*(v-4);
  516.     cx=cos(t)*r;
  517.     cz=sin(t)*r;
  518.     ay+=tt*(v-4);
  519.   }
  520. }
  521.  
  522. /*-----------------------------------------------------------------
  523.     These are some macros to ease 3D vector computation
  524. */
  525. #define diff3d(A,B,C) { \
  526.   (A##x) = (B##x) - (C##x);   \
  527.   (A##y) = (B##y) - (C##y);   \
  528.   (A##z) = (B##z) - (C##z);   \
  529. }
  530. #define cross3d(A,B,C) { \
  531.   (A##x) = (B##y)*(C##z) - (C##y)*(B##z);   \
  532.   (A##y) = (B##z)*(C##x) - (C##z)*(B##x);   \
  533.   (A##z) = (B##x)*(C##y) - (C##x)*(B##y);   \
  534. }
  535. #define len3d(A) hypot(A##x,hypot(A##y,A##z))
  536. /*-----------------------------------------------------------------
  537.    compute brightness of a triangle using light source
  538.    return a positive number less than range/2
  539. */
  540. int light_source(double px, double py, double pz,point3d *save,int range)
  541. {
  542.   static double ax,ay,az,bx,by,bz,cx,cy,cz,t;
  543.   /* get two vectors: a,b in the plane, and compute cross product */
  544.   diff3d(a,p,save[0].);
  545.   diff3d(b,p,save[1].);
  546.   cross3d(c,a,b);          /* c = a x b, normal to the plane */
  547.   t = len3d(c);
  548.   if (t<EPS) cy=1;
  549.   else { cx/=t; cy/=t; cz/=t; }  /* make unit vector */
  550.  
  551.   t = cx*lightx + cy*lighty + cz*lightz;     /* dot product = cos theta */
  552.   if (polar) t=-t;   // why?
  553.   if (range>16) range/=2;     // if large spectrum, use half
  554.   if (t>0) return range*((1-cont)*t+cont);
  555.   return range*(cont*t+cont);
  556. }
  557.  
  558. #define rline3(a,b,c,u,v,w) \
  559.   view3d(a,b,c,&x,&y);  \
  560.   view3d(u,v,w,&z,&yy); \
  561.   rline(x,y,z,yy)
  562. /*------------------------------------------------------------------
  563.    3D Graphics
  564.  
  565.    This draws a 3d surface function where the altitude is z = fz,
  566.    and x is the first variable found in fx and y is fy (if specified)
  567.    or if fy is NULL then y is the second variable found in fz.  If
  568.    there are insufficient variables in fz, then a zero is returned.
  569.  
  570.    Implementation: z is fz, x is ivar[0],y is fy or ivar[1].
  571.  
  572.    P.S.  I changed my mind in the middle of doing this.  Now y is fz,
  573.    a function of x and z, and fy is unused.  I did this because y is the
  574.    vertical axis.  Thus the code is confusing at places.
  575.    Actually, I decided to make x=fy, so that I can make a bagel in polar
  576.    mode.  So the code is very confusing at places.
  577. */
  578. int graph3d(node *fy,node *fz) {
  579.   int c,i,j,mm=0;        // mm 0=norm, 1=camera, 2=grid
  580.   static double x,y,z,xx,yy,cc,zz,x2,y2;
  581.   static point2d *save2 = 0;
  582.   static point3d *save3 = 0;
  583.   static struct { int x1,y1,x2,y2,x3,y3,x4,y4; } f;  // used by fillpoly
  584.   static char ss[60];
  585.   static int np=0,onp=0;
  586.   static double zmag,uo,us;
  587.   static short fillpat[] = { 10,11,3,8,7,900 };
  588.   int spd=0, palrot=0, palix=0, cmn,cmx;
  589.   double ca=1;
  590.   int pp=0, palst=pst, palsz=psz;
  591.  
  592.   while (1) {
  593.     if (spd) {    // speedy mode (cleardevice is very slow)
  594.       delay(50);
  595.       setcolor(0);
  596.       rline3(-2,0,0,2,0,0);         // erase axes (draw black)
  597.       rline3(0,-2*zmag,0,0,2*zmag,0);
  598.       rline3(0,0,-2,0,0,2);
  599.       rline3(2.2,0,0,2,0,-0.1);   // arrows
  600.       rline3(2.2,0,0,2,0, 0.1);
  601.       rline3(0,0,2.2,-0.1,0,2);
  602.       rline3(0,0,2.2, 0.1,0,2);
  603.       rline3(0,2.2*zmag,0,-0.1,2*zmag,0);
  604.       rline3(0,2.2*zmag,0, 0.1,2*zmag,0);
  605.     }
  606. //    if (polar) { uo=-M_PI; us=-2*uo; to=uo/2; ts=us/2; }
  607.     if (!np) {                          // reset defaults
  608.       zmag=1; us=ts; uo=to; np=25;
  609.       ca=1.0; palrot=palix=0;
  610.     }
  611.     computematrix();
  612.     xz = getmaxx()/xs;
  613.     yz = getmaxy()/ys;
  614.     if (np<3) np = 3;
  615.     if (np>2000) np = 2000;
  616.  
  617.     if (!spd) cleardevice();
  618.     spd=0;
  619. //    setlinestyle(DOTTED_LINE,0,1);
  620.     setcolor(co2);
  621.     rline3(-2,0,0,2,0,0);         // draw axes
  622.     rline3(0,-2*zmag,0,0,2*zmag,0);
  623.     rline3(0,0,-2,0,0,2);
  624.     rline3(2.2,0,0,2,0,-0.1);       // arrows
  625.     rline3(2.2,0,0,2,0, 0.1);
  626.     rline3(0,0,2.2,-0.1,0,2);
  627.     rline3(0,0,2.2, 0.1,0,2);
  628.     rline3(0,2.2*zmag,0,-0.1,2*zmag,0);
  629.     rline3(0,2.2*zmag,0, 0.1,2*zmag,0);
  630.  
  631.     if (!kbhit()) {
  632.       setcolor(co1);
  633.       sprintf(ss,"%s %g  ",msg[18],step );
  634.       if (polar==1) strcat(ss,msg[19]);
  635.       if (polar==2) strcat(ss,msg[37]);
  636.       strcat(ss,msg[36]);
  637.       outtextxy(0,0,ss);    i=0;
  638.       sprintf(ss,"%s =[%8.4f,%8.4f]  ",ivar[i].name,to,to+ts);
  639.       outtextxy(0,++i*10,ss);
  640.       sprintf(ss,"%s =[%8.4f,%8.4f]  ",ivar[i].name,uo,uo+us);
  641.       outtextxy(0,++i*10,ss);
  642.       while (i<nv) {
  643.         sprintf(ss,"%s = %8.4f  ",ivar[i].name,ivar[i].value);
  644.         outtextxy(0,++i*10,ss);
  645.       }
  646.       x = 180/M_PI;
  647.       sprintf(ss," %8.4f %8.4f %d  ",zmag,focalen,np); outtextxy(0,++i*10,ss);
  648.       sprintf(ss," %8.4f %8.4f %8.4f  ",cx,cy,cz);  outtextxy(0,++i*10,ss);
  649.       sprintf(ss," %8.4f %8.4f %8.4f  ",ax*x,ay*x,az*x); outtextxy(0,++i*10,ss);
  650.       if (mm) outtextxy(0,++i*10,msg[27]);
  651.     } else spd=1;
  652.     if (!kbhit()) {
  653.       setcolor(co2);
  654.       rline3(2,0,0.1,2,0,-0.1);
  655.       rline3(0.1,2*zmag,0,-0.1,2*zmag,0);
  656.       rline3(0.1,0,2,-0.1,0,2);
  657.       rline3(1,0,0.1,1,0,-0.1);                  // axis ticks
  658.       rline3(-1,0,0.1,-1,0,-0.1);
  659.       rline3(0.1,zmag,0,-0.1,zmag,0);
  660.       rline3(0.1,-1*zmag,0,-0.1,-1*zmag,0);
  661.       rline3(0.1,0,1,-0.1,0,1);
  662.       rline3(0.1,0,-1,-0.1,0,-1);
  663.     }
  664.     if (!kbhit() && !polar) {
  665.       setcolor(co3);
  666.       rline3(to   ,0,uo   ,to+ts,0,uo   );
  667.       rline3(to+ts,0,uo   ,to+ts,0,uo+us);
  668.       rline3(to+ts,0,uo+us,to   ,0,uo+us);
  669.       rline3(to   ,0,uo+us,to   ,0,uo   );
  670.  
  671.       rline3(-1,-zmag, 1, 1,-zmag, 1);    // draw a unit box
  672.       rline3( 1,-zmag, 1, 1, zmag, 1);
  673.       rline3( 1, zmag, 1,-1, zmag, 1);
  674.       rline3(-1, zmag, 1,-1,-zmag, 1);
  675.       rline3(-1,-zmag,-1,-1,-zmag, 1);
  676.       rline3( 1,-zmag,-1, 1,-zmag, 1);
  677.       rline3(-1, zmag,-1,-1, zmag, 1);
  678.       rline3( 1, zmag,-1, 1, zmag, 1);
  679.       rline3(-1,-zmag,-1, 1,-zmag,-1);
  680.       rline3( 1,-zmag,-1, 1, zmag,-1);
  681.       rline3( 1, zmag,-1,-1, zmag,-1);
  682.       rline3(-1, zmag,-1,-1,-zmag,-1);
  683.     }
  684. //    setlinestyle(SOLID_LINE,0,1);
  685.     setcolor(co1);
  686.     if (onp<np || onp>np*3) { onp = np;
  687.       if (save2) free(save2);
  688.       if (save3) free(save3);
  689.       save2 = malloc((np+1)*sizeof*save2);
  690.       save3 = NULL;
  691.     }
  692.     if (rain==2 && !save3) save3 = malloc((np+1)*sizeof*save3);
  693.     /*-----------------------
  694.          MAIN DRAW LOOP
  695.     */
  696.     srand(2); cmn=32767; cmx=-cmn;    // set color min/max
  697.     if (!kbhit())
  698.     for (j=0,x=to; j<=np; ++j,x+=ts/np) {
  699.       ivar[0].value = x;
  700.       for (i=0,y=uo; i<=np; ++i,y+=us/np) {
  701.         ivar[1].value = y;
  702.         xx = eval(fy);       /* by default fy = ivar[0] */
  703.         zz = y;
  704.         yy = cc = eval(fz);
  705.         if (polar==1) {       /* spherical */
  706.           x2 = xx;
  707.           xx = cos(y)*cos(x2)*yy;
  708.           zz = sin(y)*cos(x2)*yy;
  709.           yy *= sin(x2);
  710.         }
  711.         else if (polar==2) {     /* cylindrical */
  712.           x2 = xx;
  713.           xx = yy*cos(y);
  714.           zz = yy*sin(y);
  715.           yy = x2;
  716.         }
  717.         yy *= zmag;
  718.         view3d(xx,yy,zz,&x2,&y2);
  719.         x2 = xz*(x2-xo);
  720.         y2 = yz*(ys+yo-y2);
  721.         c = ca*zmag*cc*np/ts;    /* color or fill pattern */
  722.         if (c<cmn) cmn=c;
  723.         if (c>cmx) cmx=c;
  724.         if (rain==1) {          /* color by height */
  725.           c = pmod(c,palsz)+palst;
  726.           setcolor(c);
  727.         }
  728.         else if (rain==2 && i && j) {  /* color by derivative */
  729.           c = light_source(xx,yy,zz,save3+i-1,palsz) + palst;
  730.           setcolor(c);
  731.         }
  732.         if (fill) {
  733.           f.x3=save2[i].x2; f.y3=save2[i].y2;
  734.           f.x2=x2;  f.y2=y2;
  735.           if (i && j) {
  736.             if (rain && fill>2) setcolor(7);
  737.             if (rain) setfillstyle(1,c);
  738.             else setfillstyle(fillpat[pmod(c,6)],co1);
  739.             if (fill&1) fillpoly(4,&f.x1);
  740.             else fillpoly(3,&f.x2);
  741.           }
  742.           f.x1=f.x2; f.y1=f.y2;
  743.           f.x4=f.x3; f.y4=f.y3;
  744.         }
  745.         else {
  746.           if (i) line(save2[i-1].x2,save2[i-1].y2,x2,y2);
  747.           if (j) line(save2[i].x2,save2[i].y2,x2,y2);
  748.         }
  749.         save2[i].x2 = x2;
  750.         save2[i].y2 = y2;
  751.         if (rain==2 && save3) {
  752.           save3[i].x = xx;
  753.           save3[i].y = yy;
  754.           save3[i].z = zz;
  755.         }
  756.         if (kbhit()) break;
  757.       }
  758.       if (kbhit()) break;
  759.     }
  760.     /*  Draw a color legend at the bottom of the screen. */
  761.     if (rain && palsz>16 && !kbhit()) {
  762.       c = getmaxy();
  763.       for (i=0; i<256; ++i) { putpixel(i,c,i); putpixel(i,c-1,i); }
  764.     }
  765. next_key:
  766.     if (rain && palrot && v256) while (!kbhit()) {
  767.       palix += palrot;
  768.       if (palst!=16 || palsz!=240) wave_palette();
  769.       palst=16; palsz=240;
  770.       c = clock();
  771.       set_palette(palix);
  772.       while (c <= clock() && c+2 > clock());
  773.     }
  774.     c = getch(); if (!c) c=getch();  // get ONE keystroke
  775.     while(kbhit()) c=getch();        // flush kbd buffer
  776.     if (c==27) return 1;    /* escape */
  777.     if (c=='g') return 0;   /* 2d graphics */
  778.     if (c=='d') np=0;       // force reset above
  779.     if (c<84 && c>70) c+=mm;
  780.     switch (c) {
  781.     case 'c': rain=(rain+1)%3; break;     // color mode
  782.     case 'x':                         // palette mode
  783.       if (!v256) break;
  784.       pp = 1-pp;
  785.       if (pp) wave_palette();
  786.       else gray_palette();
  787.       set_palette(palix);
  788.       if (palst==16 && palsz==240 && j>np) goto next_key;
  789.       palst=16; palsz=240;
  790.       break;
  791.     case '=': palrot -= (palrot>-10); goto next_key;
  792.     case '-': palrot += (palrot<10); goto next_key;
  793.     case 'z': if (cmx==cmn) ca=1.0; else ca*=1.0*palsz/(cmx-cmn); break;
  794.     case 'f': fill=(fill+1)%(5-!rain*2); break;
  795.     case 13: mm=200-mm; break;
  796.     case 75: spincam(3); break;       // l,r,u,d
  797.     case 77: spincam(5); break;
  798.     case 72: spincam(1); break;
  799.     case 80: spincam(-1); break;
  800.     case ',': focalen*=rr;
  801.     case 82: cx*=rr; cy*=rr; cz*=rr; break;    // ins,del
  802.     case '.': focalen/=rr;
  803.     case 83: cx/=rr; cy/=rr; cz/=rr; break;
  804.     case 71: np/=rr; break;              //  home end
  805.     case 79: if (np>3) np*=rr; break;
  806.     case 73: zmag/=rr; break;   //  pgup dn
  807.     case 81: zmag*=rr; break;
  808.     case 275: ay+=tt; break;       // mm l,r,u,d
  809.     case 277: ay-=tt; break;
  810.     case 272: ax+=tt; break;
  811.     case 280: ax-=tt; break;
  812.     case 282: movecam(3); break;   // mm ins,del
  813.     case 283: movecam(-3); break;
  814.     case 273: az+=tt; break;      // mm pgup dn
  815.     case 281: az-=tt; break;
  816.     case 115: to-=ts/4; break;   // ctrl l r u d
  817.     case 116: to+=ts/4; break;
  818.     case 141: uo+=us/4; break;
  819.     case 145: uo-=us/4; break;
  820.     case 146: ts/=rr; us/=rr;
  821.               to-=ts*rr1; uo-=us*rr1; break; // ctrl ins del
  822.     case 147: to+=ts*rr1; uo+=us*rr1;
  823.               ts*=rr; us*=rr; break;
  824.     case 132: us/=rr; uo-=us*rr1; break;   // ctrl pgup dn
  825.     case 118: uo+=us*rr1; us*=rr; break;
  826.     case '?': cont-=.1; if (cont<0) cont=0; break;  // f5
  827.     case '@': cont+=.1; if (cont>1) cont=1; break;  // f6
  828.     case 'T': c='t';
  829.     default: handlekey(c);
  830.     if (c=='a')
  831.       if (polar) {
  832.         to=uo=-M_PI; ts=us=-2*uo;
  833.         if (polar==1) { to=uo/2; ts=us/2; }
  834.       }
  835.     }
  836.   }
  837.   /* never here */
  838. }
  839.